On the eigenmodes of compact hyperbolic S-manifolds 
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We present a simple algorithm for finding eigenmodes of the Laplacian for arbitrary compact 
hyperbolic 3-manifolds. We apply our algorithm to a sample of twelve manifolds and generate a list 
of the lowest eigenvalues. We also display a selection of eigenmodes taken from the Weeks space. 
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Eigenmodes of the Laplace operator contain a wealth of 
information about the geometry and topology of a man- 
ifold. This is especially true for hyperbolic 3-manifolds, 
where the Mostow-Prasad rigidity theorem [|l| ensures 
that distinct manifolds have distinct eigenvalue spectra: 
To echo the words of Marc Kac [g, one can "hear the 
shape" of a hyperbolic drum. 

Unfortunately, the eigenmodes of a compact hyper- 
bolic 3-manifold (CHM) cannot be expressed in closed 
analytic form, so numerical solutions must be sought. A 
variety of numerical methods exist to solve the problem, 
including variational principles based on the finite ele- 
ment method [|[. Perhaps the most powerful method is 
the boundary element method developed by Aurich and 
Steiner Q. Here we present an alternative method that 
came out of our work on multi-connected cosmological 
models js],^. While technically inferior to the boundary 
element method, our approach is better suited to study- 
ing a large sample of manifolds since the only inputs are 
the group generators. In contrast, the boundary element 
method requires some human effort prior to each numer- 
ical evaluation. To-date, the boundary element method 
has only be applied to two 3-dimensional examples, a 
tetrahedral orbifold 0, and the Thurston manifold [||. 
Consequently, the majority of the eigenvalue spectra de- 
scribed in this paper are completely new. 



I. PRELIMINARIES 



We seek to solve the eigenvalue problem 



A^q = qHq, 



(1.1) 



for compact hyperbolic 3-manifolds S = H^ /T, where 
the fundamental group, F, is a discrete subgroup of 
5*0 (3, 1) = PSL{2^ C) acting freely and discontinuously. 
The metric on the universal cover H^ can be written in 
spherical coordinates (p, ^, (j)): 



ds^ = dp^ + sinh^ p {dO'^ + sin^ 
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In these coordinates the Laplace operator acting on a 
scalar function ^ takes the form 
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In the simply connected space R^ ^ the eigenvalues take 
all values in the range c^ = [l,oo), and the eigenmodes 
are given by 



Qgem{p,0,<P) = X'^{p)Ye^{e,<P). 



(1.4) 



Here the Yarn's are spherical harmonics and the radial 
eigenfunctions are given by the hyperspherical Bessel 
functions 



X'Jp) = 
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The wavenumber, /c, is related to the eigenvalues of the 
Laplacian by 
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The modes have wavelength A = 27r/k and an amplitude 
that decays as l/sinh(p). The eigenmodes satisfy the 
delta-function normalization 
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In principle, the eigenmodes in the multi-connected, com- 
pact space S can be lifted to the universal cover and 
expressed in terms of the eigenmodes of H^: 
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The modes ^ q must satisfy the property 

^g(x) = "^qi^gx) \/ g eV and M x e H^ , 



(1.8) 



(1.9) 



which places restrictions on the expansion coefficients 
aqirn- Indeed, it will only be possible to satisfy (1.9) 
when q^ is an eigenvalue of the compact space. To find 
the eigenmodes we numerically solve (jL^) using a singu- 
lar value decomposition. 



II. THE NUMERICAL METHOD 



Our approach to solving (|L9|) is completely straight- 
forward. We begin by randomly selecting a collection of 



d points inside the Dirichlet domain. Taking the face- 
pairing generators ga of F, we find all images of our col- 
lection of points out to some distance Pmax in the cover- 
ing space. How we chose this distance will be explained 
later. Each point pj yields Uj images and nj{nj + l)/2 
equations of the form 



Pmax is dictated by the structure of the radial eigenfunc- 
tions Xl{p). In broad outline, they are of the form 



xiip) 



* i9o.Pj ) - * (gm ) = 0, (a ^13). 
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At each point x = g^Pj the function ^{x) is decomposed 
into eigenmodes of the covering space with wavenumber 



^fe(^) = X] X] CLMmQumix) . 



(2.2) 
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How we choose L will be explained in due course. Find- 
ing the Qkim^^ at each point is easy as there exist nu- 
merically stable recursion relations for both the hyper- 
spherical Bessel functions and the spherical harmonics. 
Schematically, we arrive at the system of equations 



The constants ^o and po depend on k and I. For fixed /c, 
po increases monotonically with increasing L Therefore, 
if we restrict our attention to some finite region with p < 
po, we need only consider a finite number of multipoles 
L. With £ held fixed, po decreases monotonically with 
increasing k. Thus, if we hold Pmax ^ Po fixed, we must 
increase L as k increases. Alternatively, if we hold L 
fixed, we must decrease pmax as k increases. 

Because the number of computational steps scales as 
L^, it makes sense to keep L as small as possible. How- 
ever, small values of L yield small values of po, which 
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The number of rows, M, is equal to 

d 

M = ^n,(n,+l)/2, 
and the number of columns, N, is equal to 

L L 
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in turn limits the number of images we are able to collect 
for each point. As a compromise, we choose L to be as 
small as possible, while leaving Pmax large enough for 
there to be at least 10 images of each point. Numerically 
we set Pmax to be the first solution to the transcendental 
equation 



Xt{p)smh{p)=0.25. 



(2.9) 
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To fix the akim^s up to an overall normalization requires 
M to equal N — 1. When M = A^ — 1 , a solution exists for 
any k. However, if M > A^ a solution only exists when k 
corresponds to an eigenvalue of the compact space. The 
standard numerical method for solving over-constrained 
systems of equations is the Singular Value Decomposition 
(SVD) H]. For a system of the form 



We also found it advantageous to fix a minimum radius 
in the same way, but with L replaced by ^min in the above 
equation. The inner cut-off helps to keep the Qum^^ of 
similar size. Finally, we found the optimal degree of over- 
sampling to lie in the range 10 -^ 100, with the high end 
of the range being required at low k and L. At higher k 
the radial eigenfunctions look less alike and a lower over- 
sampling can be used. Good all purpose choices covering 
the range /c = 1 to /c = 20 are: 



A-a = 0, 
the SVD returns solution vectors a that minimize 
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(2.7) 



L = 10 + [fc] , 
c = 10 + [100/fc] , 



(2.10) 



By incrementing in /c, the eigenvalues are revealed by 
minima in the function x^(^)- Eigenmodes with mul- 
tiplicity greater than unity will yield multiple solution 
vectors a. 

All that remains to be done is to decide on optimal 
choices for L = ^max, the tiling radius Pmax and the de- 
gree of over-constraint c = M/N. The choice of L and 



though the algorithm performs well for a wide range of 
inputs. The eigenmodes derived using the above choice 
of inputs are good to within a few percent. Choosing L 
and c larger improves the eigenmodes, but at the cost of 
slowing down the computations. 

The SVD returns solution vectors for the ak£m normal- 
ized such that 



^ ^ \akem\^ = 1- 



(2.11) 
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Thus, the eigenmodes are automaticahy delta-function 
normahzed in H^. To normahze the modes in the com- 
pact space S we numericahy perform the integral 



<i>k{x)n'{x)dV. 



(2.12) 



When k = k^ we normalize to one, and when k ^ k' we 
check that the integral vanishes (or is at least tolerably 
small). 



A. Example 
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FIG. 1. The x^ spectrum for ?7il88( — 1, 1) in the range 
A:= 1 ^ 10. 



We begin by choosing an example from Jeff Weeks' 
SnapPea [fol census of closed hyperbolic 3-manifolds and 
ask SnapPea for the face-pairing generators. Taking 
7Til88(— 1, 1) as a randomly chosen example, we ran our 
code out to /c = 10 to produce the x^(^) shown in Fig. 1. 
The x^'s for the five best solution vectors are shown. The 
eigenvalues appear as clear minima in the x^ of the best 
fit solution vector. The first mode to have multiplicity 
greater than one occurs at k = 8.34. However, below 
k = 8 we see that the x^ of the second best solution vec- 
tor occasionally exhibits minima at local maxima of the 
best fit solution vector. This occurs when the second best 
fit is formed from a weighted superposition of adjacent 
eigenmodes. In contrast, when an eigenmode has multi- 
plicity greater than one the minima of the second best 
fit coincides with that of the best fit. This behaviour is 
due to the SVD returning an orthonormal set of solution 
vectors. 

When successive minima are close together, such as 
occurs near k = 9.8 in Fig. 1, care must be taken in 
locating successive minima as they can be displaced from 
their true positions. This problem is a familiar one for 
astronomers who study spectral lines in starlight, and 
we were able to use the same deconvolution techniques 
to study our eigenspectra. The results of our analysis 
are compiled in Table I, where the eigenvalues and their 
multiplicities are recorded. 



TABLE I 


. Eigenvalue spect 


rum, q^ 


for m 


188( 


-1,1). 


20.4 

1 


22.6 

1 


27.2 

1 


30.2 
1 


39.6 
1 




46.2 
1 


51.8 
1 


55.3 

1 


60.1 
1 


70.6 
2 


75.5 
2 




78.8 
1 


80.9 
1 


83.1 
1 


86.0 
1 


96.8 
2 


98.0 
1 




99.4 
1 



III. RESULTS AND CHECKS 

We applied a battery of tests to our results. The first 
was to reproduce Inoue's [^ results for the Thurston 
manifold. Calling up 772003(— 2, 3) from the SnapPea cen- 
sus, we generated the collection of low-lying eigenmodes 
recorded in Table 11. 

Inoue was able to generate modes with q^ < 100, and 
in this range our eigenvalues agree. However, we dis- 
covered that Inoue had missed eigenmodes at q'^ = 46.2 
and q'^ = 59.1. Past q'^ = 100 we are in uncharted ter- 
ritory, so other checks have to be applied. One simple 
check is to compare the spectral staircase, i.e. the num- 
ber of modes with wavenumber < /c, with the prediction 
of Weyl's asymptotic formula 
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TABLE II. 


Eigenvalue spectrum for 


the Thurston space. 


29.3 
1 


33.5 

1 


46.2 
2 


47.8 
1 


50.8 
1 


59.1 
2 


68.9 
1 


73.8 
1 


76.2 
1 


85.8 
2 


95.1 
1 


98.0 
1 


100.1 
1 


107.5 
1 


113.8 
1 


115.5 
2 


117.4 
1 


123.3 
2 


130.3 
1 


137.9 
2 


140.0 
1 


144.5 
2 


149.8 
1 


156.3 
1 


160.5 
1 


164.1 
2 


166.6 
2 


169.5 
1 


175.0 
1 


178.2 
1 


184.6 
2 


192.3 
2 


197.8 
4 


204.3 
2 


207.2 
1 


209.8 
1 



From Fig. 2 we see that the spectral staircase closely 
follows the Weyl formula. Another check we a pplie d to 
the eigenmodes was to evaluate the integral ( |2.12D for 
modes with unequal k. The typical overlap was found to 
be on the order of a few percent, which is comparable to 
the performance of the boundary element method. 




FIG. 2. The spectra staircase for m003(— 2,3) showing its 
convergence toward the Weyl asymptotic formula. 



A. The Weeks Space 

To really put our code to the test we decided to study 
the Weeks space. This space has the distinction of being 
the smallest known hyperbolic 3-manifold, and is now 
generally thought to be the smallest example. The Weeks 
space provides a challenge to eigenmode solvers as it has 
an unusually large symmetry group. Consequently, many 
of the modes will be highly degenerate. The symmetry 
group is the Dihedral group of order 6, which has the 
presentation 



Dq =< a^b 



|a^ 6^ ba-^ba-^ > 



(3.2) 



The geometrical interpretation is that there exists a 
closed geodesic about which the space has a six-fold ro- 
tational symmetry, and a second closed geodesic, orthog- 
onal to the first, about which the manifold has a re- 
flection symmetry. If we were to choose the basepoint 
of our Dirichlet domain at one of the points where the 
two geodesies intersect, the resulting fundamental poly- 
hedron would enjoy the full Dq symmetryjj. Because of 
the Dq symmetry, we expect to find modes with 1, 2, 3, 
4 and 6 fold degeneracy. 



Taking the Weeks space m003(— 3, 1) from the SnapPea 
census, we located the first 74 eigenmodes. These are 
listed in Table III along with their multiplicities. Many of 
the higher eigenmodes are highly degenerate, in keeping 
with our expectations. The spectral staircase shown in 
Fig. 3 is in excellent agreement with Weyl's asymptotic 
formula. 

For those curious to see what the eigenmodes them- 
selves look like, we display a series of slices through a 
selection of modes. These appear in Figures 5 and 6. For 
reference we also display a view of the Dirichlet domain 
in Fig. 4, to help make contact with the 3-dimensional 
structure of the modes. 




FIG. 3. The spectra staircase for ?7i003(— 3, 1) showing its 
convergence toward the Weyl asymptotic formula. 



TABLE III. Eigenvalue spectrum for the Weeks space. 



*We thank Jeff Weeks and Craig Hodgson for providing us 
with this description. 



27.8 
1 


32.9 
2 


43.0 
1 


59.7 
1 


66.3 
1 


67.6 
2 


69.7 
2 


84.4 
1 


90.5 
1 


93.9 
1 


97.8 
2 


106.9 
2 


109.4 
2 


116.6 
1 


118.3 
1 


127.3 

1 


132.8 
1 


137.7 
4 


145.2 
2 


149.6 
2 


160.0 
1 


163.8 
1 


175.5 
3 


186.0 
2 


190.9 
2 


192.6 
3 


194.1 
2 


209.5 
1 


215.3 
4 


221.8 
2 


226.9 
2 


229.6 
3 


241.6 
1 


247.6 
1 


250.3 
1 


253.6 
2 


256.4 
6 


261.2 
1 


264.3 

1 


268.8 
1 


276.2 
1 


280.6 
2 




FIG. 4. A Dirichlet domain for the Weeks space shown in 
Klein coordinates. Like-colored sides of the polyhedron are 
topologically identified. 




FIG. 5. The first three eigenmodes of the Weeks space. The 
three views in each panel are, respectively, the x = 0, y = 
and z = slices through the fundamental domain. Here we 
are using Poincare coordinates. 




FIG. 6. The first of the 3 and 6 fold degenerate modes of 
the Weeks space. Also shown is the highest mode we gener- 
ated. 



B. The GOE Prediction 

Compact hyperbolic manifolds provide the archetypal 
setting for chaos. Consequently, we expect the statistical 
properties of the modes to be described by random ma- 
trix theory [0. Because the modes are associated with 
time-reversal invariant dynamics, we expect the statisti- 
cal properties to be those of a Gaussian Orthogonal En- 
semble (GOE). The GOE prediction is that the quantity 



\Ctkim — Ctk\ 



(3.3) 



should behave as a pseudo-random number with proba- 
bility distribution 



P{x) 



^27TX 



-x/2 



(3.4) 



In the above equations, ak denotes the average of the 
akim^s and a^ their variance. To avoid the singularity at 
X = 0, it is conventional to compare numerical results to 
the cumulative distribution 
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P{x) dx = eri{yx/2) . 



(3.5) 



Taking the mode with eigenvalue q^ = 175.5 as an 
example, we display in Fig. 7 the first 676 akim^s in a 



scatter plot. Notice that the distribution is independent 
of £ and m, as expected for a chaotically mixed state. 
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FIG. 7. The a/c^m's for the eigenmode q^ = 175.5, shown 
up to imax = 25. The numbering scheme for the modes maps 
dim ^-^ cii where i = i'^-\-i-\-l-\-m. 



The cumulative distribution I{x) is compared to the 
GOE prediction in Fig. 8. The agreement is quite re- 
markable. 




FIG. 8. The cumulative distribution I{x) found for the 
eigenmode with q^ = 175.5 (blue line) and the GOE pre- 
diction (red line). 



C. The circles test 

The next test we applied to the modes is one that 
we plan to apply to our own universe H,^,|l2[. Imagine 
drawing a 2-sphere of radius p about the basepoint of 
the Dirichlet domain. Viewed in the universal cover, the 
space will contain an infinite number of copies of this 



2-sphere. If the radius of the 2-sphere exceeds the in- 
radius of the Dirichlet domain, then the 2-spheres will 
intersect along a circle. Mapping the entire picture back 
inside the Dirichlet domain, we see that the 2-sphere self- 
intersects. If we now take a 2-sphere slice through one of 
the eigenmodes, the amplitude of the mode must match 
up around the matched pair of circles. 

The largest matched circles lie on face-planes of the 
Dirichlet domain. Taking a face-pairing generator ^, and 
representing it as a 4 x 4 real matrix in 5'0(3, 1), the 
angular radius a of the matched circle is 



a = arccos 



^00 - 1 



V^oo - 1 tanh(p) J 
The (6>, (j)) coordinates of the circle centers are 

f arccos ( gio/ygoo ~ ^ ) ' arctan (^30/^20) j , 



(3.6) 
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We display in Fig. 9 the amplitude of the eigenmodes 
q^ — 43.0 and c^ = 84.4 around a 2-sphere of radius p = 1 
using an equal-area projection. Two of the matched circle 
pairs are indicated by white lines. In Fig. 10 we plot the 
amplitude of each mode around each pair of circles and 
see that they are indeed properly matched. If our uni- 
verse is multi-connected, we hope to see^ similar matched 
circles in the cosmic microwave sky 




FIG. 9. A two-dimensional slice through the c^ — 43.0 (up- 
per panel) and cf' — 84.4 (lower panel) eigenmodes. The slice 
is taken on a 2-sphere of unit radius. Two pairs of matched 
circles are also indicated. The relative phasing of one matched 
pair is shown on the lower panel. 
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FIG. 10. Each panel shows ^q around a pair of matched 
circles, C (in red) and C (in blue). Also shown (in green) is 
the difference ^g — ^q around the matched pair. This small 
difference represent the numerical error in our algorithm. The 
first panel correspond to the mode q^ — 43.0 sampled along 
the inner pair of circles. The second panel shows the same 
pair of circles, but for the mode q^ = 84.4. The third and 
fourth panels are likewise shown for the outer pair of circles. 



IV. THE LOWEST EIGENMODES 

Having established that our algorithm is reliable, we 
set it loose on the SnapPea census to produce the list 
of lowest eigenvalues recorded in Table IV. The lowest 
eigenvalue, g^, is a useful topological invariant that has 
attracted considerable attention in the mathematical lit- 
erature. There exist a variety of upper and lower bounds 
on qf. A summary of these bounds can be found in the 
works of Callahan [|l3[ and Cornish et al. [|lj]. Amongst 
the sharpest are those that employ the diameter^], D^ of 
the space: 



AD 



D^{smh{D)^Dy 



< 



Qi 



<1+ ^ 
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Here D denotes the square root of the smallest integer 
that is greater than or equal to D^ . A listing of the di- 
ameters can be found in Table V. In all cases, the lowest 
eigenvalu es fo und by our algorithm fell in the range dic- 
tated by (4T). The eigenvalue bound ( |4j] ) tells us that 
the wavelength of the lowest eigenmode must be greater 
than or equal to the diameter of the space. Curiously, 
we found that Ai = (1.3 ^ 1.6)1) for all twelve examples 
studied. 

On a cautionary note, the eigenvalues listed in Table 
IV might not be the lowest supported by these spaces. 
Our method for finding the eigenmodes is unable to de- 
tect modes with g'^ < 1, as these modes have imaginary 
wavenumbers. Moreover, the spherical eigenmodes we 
use as our expansion basis start to look very much alike 
for wavenumbers k < 1, so we may have missed modes 
in the range k = [0, 1/4]. We are currently developing a 
variant of our method to handle all modes below q^ = 2. 
At the other end of the spectrum, the only limitation in 
going out to g^ = oo is computer power. To get all the 

TABLE IV. Lowest eigenvalues 



Vol 



G 



q! 



mi 



m003(-3,l) 
m003(-2,3) 
s556(-l,l) 
m006(-l,2) 
ml88(-l,l) 
v2030(l,l) 
m015(4,l) 
s718(l,l) 
ml20(-6,l) 
s654(-3,l) 
v2833(2,3) 
v3509(4,3) 



0.9427 
0.9814 
1.0156 
1.2637 
1.2845 
1.3956 
1.4124 
2.2726 
3.1411 
4.0855 
5.0629 
6.2392 



Dg 

D2 
Z4 
D4 
D2 
D2 
D2 
D2 
Z2 
D2 
Z2 
D2 



27.8 
29.3 
27.9 
21.1 
20.4 
16.2 
28.1 
10.1 
7.50 
5.88 
6.29 
6.06 



'The diameter is defined to equal the greatest distance be- 
tween any two points in the manifold. 



TABLE V. 


L^, 


diameter, 


wavelen^ 


5th and wavenumber 


S 




L^ 


D 


\i/D ki 



m003(-3,l) 


0.585 


0.843 


1.44 


5.18 


m003(-2,3) 


0.578 


0.868 


1.36 


5.32 


s556(-l,l) 


0.831 


0.833 


1.45 


5.19 


m006(-l,2) 


0.575 


1.017 


1.38 


4.48 


ml88(-l,l) 


0.480 


0.995 


1.44 


4.41 


v2030(l,l) 


0.366 


1.082 


1.49 


3.90 


m015(4,l) 


0.794 


0.923 


1.31 


5.21 


s718(l,l) 


0.339 


1.439 


1.45 


3.01 


ml20(-6,l) 


0.314 


1.694 


1.45 


2.55 


s654(-3,l) 


0.312 


1.946 


1.46 


2.21 


v2833(2,3) 


0.486 


1.701 


1.60 


2.30 


v3509(4,3) 


0.346 


1.802 


1.55 


2.25 



modes out to q^ = 250 takes around 12 hours on a single 
RIOOO Silicon Graphics chip (or 10 minutes if you use a 
64 processor Origin 2000 supercomputer as we did). 
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